Slow relaxation and aging kinetics for the driven lattice gas 
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We numerically investigate the long-time behavior of the density-density auto-correlation function 
in driven lattice gases with particle exclusion and periodic boundary conditions in one, two, and three 
dimensions using precise Monte Carlo simulations. In the one-dimensional asymmetric exclusion 
process on a ring with half the lattice sites occupied, we find that correlations induce extremely slow 
relaxation to the asymptotic power law decay. We compare the crossover functions obtained from our 
simulations with various analytic results in the literature, and analyze the characteristic oscillations 
that occur in finite systems away from half-filling. As expected, in three dimensions correlations are 
weak and consequently the mean-field description is adequate. We also investigate the relaxation 
towards the nonequilibrium steady state in the two-time density-density auto-correlations, starting 
from strongly correlated initial conditions. We obtain simple aging scaling behavior in one, two, 
and three dimensions, with the expected power laws. 

PACS numbers: 64.60. an, 64.60. De 



I. INTRODUCTION 

Driven lattice gases represent perhaps the simplest 
models of interacting non-equilibrium systems that nev- 
ertheless display intriguingly complex features, and con- 
sequently have been carefully studied for many years 
now pHll. Since there exists to date no general theo- 
retical framework to adequately capture and classify the 
macroscopic properties of non-equilibrium systems, sim- 
ple paradigmatic models such as driven lattice gases have 
proven useful to gain a thorough understanding of at least 
certain models driven towards out-of-equilibrium station- 
ary states, characterized by non- vanishing macroscopic 
(probability, particle, or energy) currents. 

Let us briefly review the basic features of driven lattice 
gases with particle exclusion to provide a proper context 
for the present investigation. For the one-dimensional 
driven lattice gas with periodic boundary conditions, 
equal-time correlations have long been established to be 
trivial and governed by a simple product measure 0] . Its 
temporal evolution, in contrast, displays rich features, 
and is not easily accessible analytically. The Bethc ansatz 
method was used to calculate the spectral gap at half fill- 
ing which gives the value of the dynamic exponent . 
This approach was later generalized for arbitrary densi- 
ties [HI UM ■ The Bethe ansatz also allows for calculation 
of conditional time-dependent probabilities for the lattice 
gas [H, [Hi- Yet time-dependent correlation functions 
have only been amenable to very limited general exact 
solutions most of which are tractable only for very small 
system sizes. The correlation functions are governed by 
non-trivial scaling laws even in the stationary state 0). 

It is well-understood that the large-scale behavior of 
the driven lattice gas in one dimension is governed in 
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the coarse-grained continuum limit by the noisy Burgers 
equation (14] . and is also mappable to a surface growth 
model that belongs to the KPZ universality class [l5j |. 
The scaling form of the density-density correlation of the 
driven lattice gas has been derived and the associated 
scaling exponents calculated by means of renormalization 
group techniques [l4l - [l6j . Scaling functions have been 
computed by means of the mode-coup ling ap proximation 
[17H23| . exact numerical integration [24| - l26j . and other 
analytical approaches p7l - l30j . 

In this work, we employ extensive large-scale Monte 
Carlo simulations to confirm the exactly known scaling 
exponents in one, two, and three dimensions. We will 
also compare our numerical data to different predictions 
for the scaling functions. In one dimension in particu- 
lar, our Monte Carlo simulations yield a remarkably slow 
approach to the asymptotic scaling behavior, governed 
by an unusual form of the associated scaling function. 
Non-stationary properties of the driven lattice gas in the 
physical aging regime are considerably less understood. 
For KPZ surface growth, a renormalization group analy- 
sis demonstrated that the so-called initial slip exponent 
can be related to the dynamical exponent, and no addi- 
tional singularities associated with the initial-time sur- 
face emerge 3l|. We find that similar scaling relations 
hold for driven lattice gases in arbitrary dimensions. In 
addition, for appropriate highly correlated initial condi- 
tions (that in the growth model mapping correspond to 
a flat surface), we establish simple aging scaling in the 
non-stationary time regime (32j . This is in accord with 
studies that showed si mple ag ing as well for both the KPZ 
surface growth model [33l436j and its linear counterpart, 
the Edwards- Wilkinson model (36l - [39| . A few exact re- 
sults were found in the non-stationary regime, for small 
systems [4fj| - |43l |. 

The outline of this paper is as follows. In the follow- 
ing section we describe the models under investigation, 
comment on our Monte Carlo simulation specifics, de- 
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fine the relevant observables, and provide the dynamic 
scaling forms for the density-density (and height-height) 
correlation function. In section Mil we present our simu- 
lation results in the non-equilbrium steady state for the 
model, focusing on the behavior in one dimension. Short 
time results for densities away from half filling yield char- 
actaristic oscillations in the density-density correlation 
function. Next finite-size scaling features are discussed 
in some detail. We generate the scaling functions from 
our simulation data and compare them to the accepted 
forms. The unusually slow relaxation of the local expo- 
nent of the density-density correlation function towards 
its asymptotic universal value is carefully investigated for 
the one-dimensional case, and contrasted with the much 
more straightforward behavior seen in higher dimensions. 
We conclude this section with a discussion of the effects 
of finite rather than infinite hopping bias. Section IIVI 
contains our simulation results for non-stationary prop- 
erties of the driven lattice gas, starting from a correlated 
initial state, in one, two and three dimensions. We cast 
our results in the standard simple aging scaling form and 
compute the aging exponents by means of a simple scal- 
ing argument, which is a-posteriori confirmed by our sim- 
ulations. 



II. THE DRIVEN LATTICE GAS WITH 
EXCLUSION 

A. Model description 

We study driven lattice gases with site exclusion on 
a d-dimensional square lattice of size L d with periodic 
boundary conditions, with a total number of N parti- 
cles 0, y, H[. The exclusion constraint limits the site 
occupation number at a position x on the lattice to the 
values n(x) £ {0, 1}. The particles may hop to nearest- 
neighbor lattice sites only if these are empty. In one 
dimension, this model is referred to as the periodic asym- 
metric simple exclusion process (ASEP). We denote the 
probabilities of hopping "forward" and "backward" as p 
and q = 1 — p, respectively For p = q, no net parti- 
cle current arises; the system is then in equilibrium. For 
p = 1 and hence q = 0, only unidirectional motion (to 
the right) is possible, and the model reduces to the to- 
tally asymmetric simple exclusion process (TASEP). The 
generalization of the driven lattice gas model to higher 
dimensions d is straightforward. After a single specific 
drive direction is selected, we define the hopping proba- 
bilities (to empty sites) along this direction to be pu and 
qu = 1 — P\\, while there is no hopping bias in the d — 1 
transverse dimensions, p± = q± (symmetric exclusion 
process, SEP). 

In one dimension, this lattice gas model can be mapped 
to a surface growth problem in the KPZ universality class 
flol ]. To this end, one defines local Ising spin variables 
cr(x) = 1 — 2n(x), where n(x) is the occupation number. 
The height variable at site x is then given by the accumu- 



lated magnetization from the (arbitrary) reference point 
at site 1 up to x: 



(1) 
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Particle exclusion restricts the local slope of the surface 
height profile to ±1. This mapping allows us to directly 
relate observables such as correlation functions calculated 
for the surface growth model to those for our driven lat- 
tice gas model. 



B. Monte Carlo simulations 

We employ both standard Metropolis and continuous- 
time Monte Carlo algorithms to study the driven lat- 
tice gas. The standard Monte Carlo algorithm [44[ starts 
with all particles placed at random positions. A particle 
is randomly selected, and then one of the 2d directions 
for a potential nearest-neighbor hop is randomly cho- 
sen. A hop has a predetermined probability Ri 6 [0, 1] 
(i = l,...,2d, Ri is either pn, qu, p±, or q±). A ran- 
dom number r £ [0, 1] is picked; if r < Ri then hop i 
is performed. Particle exclusion is realized in this sim- 
ulation by setting Ri = for any hop that attempts to 
move a particle into an occupied site. One Monte Carlo 
time step has passed once N hopping attempts have been 
made. On average the algorithm thus attempts to hop 
each particle once per Monte Carlo time step. In dimen- 
sions d > 2, the hopping probabilities are asymmetric 
parallel to the drive (if not mentioned otherwise we use 
PU = 1, <j|| = 0), and symmetric in the transverse direc- 
tions p± = qj_ = 1/2. With the periodic boundary condi- 
tions, this prescription generates a (fluctuating) particle 
current along the drive direction. 

For the one-dimensional system we employ the 
continuous-time Monte Carlo algorithm [3, HE) , which 
is initialized the same way as the standard Monte Carlo 
process. The crucial feature of this algorithm is that hops 
are never rejected; they are simply performed at the cor- 
rect rates. The probability for each possible hop i from 
a given configuration C of the system is Ri. An active 
list is kept of all possible hops with Ri ^ 0, and the next 
hop is chosen by sampling the list. Time is incremented 
by At = 1/J2i Ri- For the driven lattice gas with exclu- 
sion, each site has an equal probability of being occupied, 
equal to the filling fraction or mean density p = N/L; 
e.g., for a half-filled system the probability is p = 1/2. 
This algorithm eliminates all attempted hops that would 
be rejected due to exclusion. This is realized in the simu- 
lation as follows for the ASEP. Two lists are maintained, 
L ri Li, one for the hops to the right and left with prob- 
abilities p and q, normalized by p + q = 1. A random 
number r G [0, 1] is generated. If r < p, a hop to the right 
is performed, if r > p we chose a hop to the left. Note 
that we may only proceed in this manner because the 
list lengths Nh of possible left and right hops are equal. 
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A random integer n G [1, A/ t ] is selected, and we corre- 
spondingly update the nth particle in the L r or L\ list. 
The time counter is increased by At which is calculated 
using the rates prior to the hop being performed. This 
continuous-time algorithm yields an efficiency enhance- 
ment of approximately a factor of 2, because only hops 
which can be performed arc chosen, and consequently al- 
lows us to increase the effective simulation time by that 
factor also. We note that in the following, all distances 
and times will be measured in units of the lattice spacing 
and Monte Carlo Steps (MCS), respectively. 



C. Quantities of interest 

Our starting point is the two-time density-density 
(connected) correlation function (cumulant), defined as 

S(x ,x,t ,t) = (n(x,t)n(x ,t )) — (n(x ,t))(n(x,t)), (2) 

where (■) denotes an ensemble average over (many) 
stochastic realizations. In the stationary state, a peri- 
odic system becomes translationally invariant in both 
space and time. Therefore the correlation function de- 
pends only on coordinate and time differences. Defining 
the average density p = (n(x,t)) = N/L, we then obtain 



equal-time correlation function becomes 



S{x' -x,t' -t)= (n(x,t)n(x" ,t')) - p 2 



(3) 



For the driven lattice gas (ASEP) in the stationary 
state, it is well known that each possible particle config- 
uration C occurs with equal probability P*(C), namely 
the inverse of the number N(C) = (^) of possible states 



P*(C) 



1 



N(C) 



Nl(L — N)\ 
L\ 



(4) 



This independent-particle result remarkably remains 
valid even in the presence of on-site particle exclusion 
However, if nearest-neighbor interactions were 
added, the stationary probabilities would depend on the 
particles' relative positions. 

^From the simple product measure ((4]) , we can readily 
obtain the equal-time stationary density-density correla- 
tion function S(x' — x^0,t' — t = 0). The probability 
to find a particle at a given site x is just the average 
density p = N/L, and given that site x is occupied the 
density at a different site x! ^ x is (N — 1)/(L — I), 
which yields (n(x,t)n(x' ,t)) = J^c P* {C)n(x)n(x') = 
J2 C > P*{C')N(N - l)/L{L - I). The sum now extends 
over configurations C, which encompass all possible 
states of N — 2 particles in the remaining L — 2 sites, 
whose stationary probability P*(C) is again given by 
Eq. But the summation gives us a multiplicative 

factor N(C') = 1/P*(C), which is the number of possi- 
ble configurations of the remaining N — 2 particles placed 
in the remaining L — 2 sites. Hence we arrive at the sim- 
ple result (n(x,t)n(x',t)) = N(N - 1) / L{L - 1), and the 



S(x' -x^0,0) = 



N N - 1 
L L-l 



(5) 



Taking the thermodynamic limit N, L — > oo while keep- 
ing p = N/L constant, the correlation function ([5]) van- 
ishes for x ^ x! . However, in a finite system, Eq. ([5]) 
yields a negative value 



(PL-I) 
L-l 



P(l ~ P) 
L-l ' 



(6) 



reflecting the particle anti-correlations induced by site 
exclusion. 

A more interesting quantity to study is the two-time 
(connected) auto-correlation function S(x, x, t', t), which 
in the stationary state can be written as 5(0, t) = 
{n(x, 0)n(x,t)) - p 2 . 

In a finite one-dimensional system with periodic 
boundary conditions, i.e., a ring with L sites, another 
notable quantity is the velocity of a single "tagged" par- 
ticle 



Vi 



L-N 
L-l 



(7) 



We have checked and confirmed this tagged particle ve- 
locity in our Monte Carlo simulations. Note that vt = 
I — p in the thermodynamic limit. Since our numerical 
data are obtained in finite systems, the typical return 
time t r = L/v t for a specific particle to traverse the en- 
tire ring and come back to its original site posits an upper 
limit to the argument t < t r for a meaningful analysis 
of the time-dependent auto-correlation function 5(0, t). 
We remark that the tagged particle velocity has to be 
carefully distinguished from the propagation speed of a 
collective density fluctuation Q 

v c = 1 - 2p, (8) 

which vanishes only at half-filling p = 1/2. 



D. Dynamic scaling 

In the scaling limit, the steady-state density-density 
correlations for the driven lattice gas with exclusion in 
d dimensions are described by a general homogeneous 
function of the following form [TJ 



S(x_L,x h t) = b d+A s(bx ± ,b 1+A x h b z ty 



(9) 



where b represents an arbitrary real scale factor, and the 
subscripts _L and || denote transverse and parallel direc- 
tions with respect to the drive. The exponents A and z 
are the anisotropy and (transverse) dynamical exponents, 
respectively. A renormalization group analysis [161 ] yields 
the upper critical dimension d c = 2, and the exact ex- 
ponent values z = 2 and A = (2 — d)/3 for d < d c . For 



4 



d > d c the exponents are given by mean- field theory, i.e., 
z = 2 and A = 0, with logarithmic corrections at d c = 2. 
Setting b = t~ x l z , Eq. (0 becomes 

S(x ± ,x h t) = r^y'uit-V'Sut- 1 '*^, (10) 

where Z|| = z/(l + A) is the dynamical exponent along 
the drive direction. Taking x±, x\\ — > 0, we obtain for the 
density auto-correlation function in the scaling limit: 



S(p,o,t)~t-<, c 



d + A 



rf + A 
(1 + A)z 



(11) 



In order to carefully characterize the approach to the 
asymptotic regime in our Monte Carlo simulations, we 
shall utilize a corresponding local (effective) exponent de- 
fined via 



CeffW 



Iog5(0,t)-log5(0,t-l) 
logt-log(t-l) 



(12) 



A more detailed discussion of each dimension is in order 
at this point: 

d = 1. The transverse spatial directions become ob- 
solete, and with the exact renormalization group result 
A = 1/3 we obtain zn = 3/2 and C = l/z|| = 2/3. As 
mentioned above, the scaling behavior for the equiva- 
lent surface growth model falls into the KPZ universality 
class [HI]. The central quantity in this description is the 
height-height correlation function 



C(x,t) = ([h(x,t) - h(0,0)f 



(13) 



which is related to the density-density correlation func- 
tion via IH 



Asymptotically, it takes the scaling form 
COM) = b- 2a C{bx,b z H), 



(14) 



(15) 



with the roughness exponent a (note that in the growth 
model literature, the dynamical exponent z\\ is usually la- 
beled z). According to Eqs. (|14[) and (fl5|) . we can there- 
fore identify 2 — 2a = (d + A)z\\/z = 1 in one dimension, 
or a = 1/2 = 2 - z\\ 

d = d c = 2. At the upper critical dimension, loga- 
rithmic corrections arise. The simple power law (|11[) is 
replaced with 0, d 



5(0,0, t) 



i(logt)V2 



(16) 



Accordingly, in order to capture the approach to the 
asymptotic scaling limit we need to redefine the local 
exponent in this case as 



S'(0,0,t)(log<) 1 / 2 <xt~^ f{it \ 



It turns out though that the mean-field approximation 

S M F(0,0,i) -t- 1 (18) 

provides reasonable agreement with our simuation data 
as well. 

d > 2. Above the critical dimension, the mean-field 
predictions for the scaling exponents A = and z = 2 = 
zh , C = d/2 should provide an adequate description of 
the density auto-correlation function data. 

Finally, we address the two-time density auto- 
correlation function S(0,t,s) in the non-stationary 
regime when the system relaxes from a strongly corre- 
lated initial state that considerably differs from any non- 
equilibrium steady-state configuration. As we shall ex- 
plain in more detail in Sec. IV. B below, one expects that 
the physical aging scaling regime s <S t is governed by 
the simple aging scaling form [32| 



5(0, t,s) 



- c h s (t/t 



(19) 



Indeed, we shall see that this scaling ansatz yields satis- 
factory data collapse. We remark that in the literature 
on glassy systems, an alternate so-called subaging fitting, 
5(0, t,s) = hsit/s^) with fj, < 1, has been popular [32j ]. 



III. STEADY-STATE RESULTS 

A. Recurrence oscillations in small systems 

We begin with a discussion of small-scale properties of 
the driven lattice gas on a one-dimensional ring of length 



o 

t3? 




(17) 



FIG. 1. (Color online.) Characteristic recurrence oscilla- 
tions in the density auto-correlation function 5*(0, t) for a 
one-dimensional driven lattice gas (TASEP) at low density 
p — 0.1 and L — 100 (close-up). The horizontal (red) line in- 
dicates the minimum value Smin = -0.909 x 10~ 3 (red line), 
see Eq. ([6]). The inset has been included to show the actual 
magnitude of the oscillations. The return time is T c = 125 
MCS. The data are averaged over 8,000,000 realizations. 
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FIG. 2. (Color online.) Recurrence oscillations in the driven 
lattice gas (TASEP) auto-correlation function S(0,t), again 
for p = 0.1, but a larger one-dimensional ring of length L — 
2000; here S m i n — 0.45 x 10 -4 , indicated by the dashed (red) 
line. The return time is T c = 2500 MCS, and the data are 
averaged over 8,000,000 realizations. 



FIG. 3. (Color online.) Fourier transform |S(cj)| of S(0, t) for 
one-dimensional driven lattice gases (TASEP) with L = 400, 
at various densities p = 0.15, 0.25, 0.35, and 0.5. The curves 
are vertically shifted to separate them and show the peaks. In 
each case, the data were averaged over 800,000 realizations. 



L (TASEP). When the total density p = N/L ^ 1/2, col- 
lective density fluctuations travel around the system with 
mean velocity v c = 1 — 2p. This causes characteristic os- 
cillations in the time-dependent density auto-correlation 
function S(0,t), with peaks at multiples of the return 
time T c = L/v c , as can be seen clearly in Figs. [T]and[21 
which show data obtained for L = 100 and L = 2000, re- 
spectively, both at low density p = 0.1. (We only report 
results for densities less than 1/2, since as a consequence 
of particle-hole symmetry, the results for p = 0.5 — 8 are 
equivalent to those for p = 0.5 + 5.) 

Observe the flat section in Fig. Q] at a negative value 
5 min = -0.9 x 10~ 3 in the first minimum. Invoking the 
fact that in the steady state, all configurations carry the 
same statistical weight, we can understand this lower 
bound on the density auto-correlation function as fol- 
lows: Consider a single occupied site x at time t, hence 
n(x,t) = 1. At time tl > t, assume that the particle 
has left the site x, but has not yet traversed the entire 
system and returned: this situation will give the lowest 
average possible value for the correlation function. While 
the new location of the particle previously at x can be 
inferred from the tagged particle velocity v t , Eq. (UJ), the 
other TV — 1 particles are left to fill any of the remaining 
L — 1 sites with equal likelihood, which gives a proba- 
bility (N — 1)/(Z — 1) for site x to become occupied at 
t' . Averaging over all L sites, we thus arrive at the same 
average minimum value © for the auto-correlation func- 
tion S(0, t > 0) when the original particle has traveled 
away that we derived earlier for the equal-time density 
cumulant for S(x* — x ^ 0, t = 0): After a brief decorrela- 
tion time, typical configurations at any given site at later 
time t' are statistically correlated to the reference con- 
figuration at t < t' in the same manner as different sites 



are at equal time in the stationary state. The amplitude 
of the recurrence oscillations naturally decreases with in- 
creasing system size. Note, however, that the mean neg- 
ative value 5 m in = 0.45 x 10~ 4 is still clearly visible in 
Fig. H for L = 2000. 

In order to extract the period of the recurrent density 
fluctuations, we Fourier transform the auto-correlation 
function, S(u)) = J S(0,t)e lult dt. The first peak in \S(u>)\ 
is caused by the movement of a density fluctuation at 
a frequency uj c = 2ir/T c , see Fig. [3] At low densities, 
higher harmonics are also visible. The characteristic fre- 
quency uj c vanishes as the density approaches 1 /2. We 
finally mention previous work by Adams et al. [48j and 
Cook and Zia [49[ who studied density oscillations in the 
open TASEP that are also caused by fluctuations that 
traverse the entire system. Gupta et al. studied similar 
oscillations in the displacement variance for tagged par- 
ticles [50| . Exact calculations of the spectral gap yield 
evidence of oscillations as well d, [Io|, El , which are also 
found in exact solutions for the time-dependent condi- 
tional probabilities and single particle correlation func- 
tion [Kit - 



B. Finite-size scaling 

Next we analyze the finite-size scaling properties of 
the density auto-correlation function for the driven lat- 
tice gas with exclusion in one dimension (TASEP). The 
expected finite-size scaling form is readily obtained from 
the general expression © by setting L = b~( 1+A \ and 
omitting the transverse space directions: 

S(x, t) = L^Sfs (x/L, t/L z ny (20) 
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FIG. 4. [Color online.) Finite-size scaling for the density 
auto-correlation function for L = 64 (data averaged over 
80,000 realizations), L = 256 (10000 realizations), L = 2048 
(40,000 realizations), and L = 2,097, 152 (60 realizations). 




FIG. 5. (Color online.) Double-logarithmic plot of the scaling 
function F(£), defined in Eq. ((22]), obtained from Monte Carlo 
simulations in a driven lattice gas on a ring of length L — 
2000, with x 6 [0,64] and t <E [0,64]. The data are averaged 
over 200 realizations. 



Taking x — > wc find for the auto-correlation function 

S(Q,t) = L- 1 S FS (t/L z ny (21) 

As shown in Fig. @J the finite-size scaling form (f2"Tj) pro- 
duces satisfactory data collapse with zn =3/2 for our 
simulation results for vastly different lattice sizes rang- 
ing from L = 64 to L = 2,097,152. At t/L z W ps 0.1 
one observes the expected finite- L cutoff from the scal- 
ing form (|21[) . as the return time t r is approached. We 
note that earlier work has studied the finite-size scaling 
properties of the TASEP with open boundaries in the 
maximum current phase which is equivalent to the peri- 
odic case; however much smaller system sizes of respec- 
tively i = 50 and L = 257 were utilized by Picrobon et 
al. [531 and Juhasz and Santen 1531. 



C. Long-time properties 

1. Scaling functions 

We now address the stationary height-height correla- 
tion function (|13|) in one dimension, and explicitly deter- 
mine associated scaling functions. Setting b = 1/|%L and 
recalling a = 1/2 and z\\ = 3/2, Eq. (fl~5j) becomes [lj| 

C(x,t)oc\x\F(£), £ = t/M 3/2 , (22) 

whereas ecjuivalently the matching condition b = i~ 1//z n 
leads to [25( 

C(x,t) oc t 2 / 3 g(y), y = x/t 2 ^ 3 . (23) 
Consistency then requires that 

F(0 <x e'\ 9(y) oc \y\ (24) 



for large arguments £ 3> 1 and i/ > 1. Figures [5] and 
[6] depict the scaling functions F(£) and g(y) as obtained 
from our simulations for a driven periodic lattice gas with 
exclusion on a ring of length L = 2000. The Monte Carlo 
data nicely confirm the relations ([24")) . The crossover from 
= const, for small arguments £ to the asymptotic 
behavior ~ £ 2 / 3 visible in Fig. [5] can be compared with 
the numerical solution of the mode-coupling equations 
depicted in Fig. 2 of Re f. [ 19t] and the simulation results 
shown in Fig. 2 of Ref. [54[ . The mode-coupling approx- 
imation appears to predict sharper crossover features, 
with F(£) remaining constant until £ s=s 0.5, where our 
data already indicate a noticeable curvature in the scaling 




y 



FIG. 6. (Color online.) The scaling function g(y) defined 
in Eq. (f23|) . from the same data as in Fig. [5] (L = 2000, 
x € [0,64], and t £ [0,64], averaged over 200 realizations). 
The inset represents the same data in a double-logarithmic 
plot to display the crossover to a constant as y — > 0. 
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FIG. 7. (Color online.) Logarithmic plot of the scaling func- FIG - 8 - ( Color online.) Double-logarithmic plot of the sta- 
tion f(y), defined in Eq. (25|, for L = 2000, with x £ [0, 256] tionary auto-correlation function 5(0, t) for a one-dimensional 
and t G [0, 256], with the data averaged over 400,000 realiza- driven lattice S as of len 8 tn L = 2048 with exclusion. The 
tions. The error bars indicate one standard deviation. data are a veraged over 40,000 realizations. The light green 

line represents a power law with exponent —2/3. 



function F. The simulation data in Rcf. 54 1, obtained 
for much larger systems with L = 1048576, also display 
sharper crossover features than our data for L = 2000. 

Combining Eqs. (fT4|) and (|23|) . we recover the scaling 
form (fit)]) for the density-density correlation function in 
one dimension, 



S{x,t)<xt-*l*f{y). 



(25) 



Our Monte Carlo simulation data for the scaling func- 
tion f(y), measured by averaging over a very large num- 
ber of realizations, are plotted in Fig. [7] Within the er- 
ror bars, our data are in reasonable aggreement with the 
stretched exponential function f(y) oc exp(— cy 3 ), where 
c = 0.295(5), that was computed in Rcf. [251 ] . 

2. Local exponent 

In Fig. [8] we plot the density auto-correlation func- 
tion S(0,t) as obtained from our simulations for a one- 
dimensional driven lattice gas with exclusion with L = 
2048 (sec also Fig. g]). From Eq. (HI]) with £ = l/z\\ = 
2/3, we expect the asymptotic long-time power law 
S(0,t) ~ t~ 2 / 3 , which fits the Monte Carlo data well for 
t > 100 on the double-logarithmic plot. The approach to 
the expected asymptotic power law auto-correlation de- 
cay can be carefully probed by studying the local effective 
exponent ( c s(t) defined in Eq. (|T2l . 

Figure [9] shows that the local exponent increases with 
time and in fact surpasses the asymptotic value 2/3, indi- 
cated by the horizontal line. In order to demonstrate that 
this is actually a finite-size effect, albeit a quite unusual 
one, we display in Fig.[10]a comparison of simulation data 
for the time dependence of the effective exponent for sys- 
tems of sizes L = 2048 and L = 1024 x 2048 = 2097152 
on a logarithmic time scale. Even in the much larger 



lattice, Ceff(i) still rises beyond the asymptotic exponent 
value £ = 2/3 of the infinite system, yet at later time. 
Since we have not observed this overshooting in unbiased 
lattice gases with exclusion, it is clearly a consequence of 
the non-equilibrium drive. We note that the extremely 
slow crossover towards the asymptotic temporal power 
law has been recorded before in a model of reconstitut- 
ing dimers that can be mapped to the TASEP [55[. 

We next compare the behavior of the local auto- 
correlation exponent Ceff(i) for the one-dimensional 
driven lattice gas with its two- and three-dimensional 
counterparts, see Figs. [11] and Q21 In each case, the 
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FIG. 9. (Color online.) The local auto-correlation exponent 
(|12[) obtained from Monte Carlo simulation data for a one- 
dimensional driven lattice gases with length L — 2048 and for 
p — 1/2. The data are averaged over 40,000 realizations. The 
error bars represent one standard deviation, and the (red) 
horizontal line the asympotic value 2/3. 
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FIG. 10. (Color online.) Time dependence of the local ex- 
ponent Ccff(^) f° r one-dimensional driven lattice gases with 
L = 2048 (black circles, data averaged over 40,000 realiza- 
tions) and L — 2097152 (red squares, 60 realizations). 



FIG. 12. (Color online.) The local exponent Ceff(i) for a 
three-dimensional driven lattice gas with L x — L y — L z — 50 
and p = 1/2. The data are averaged over 1,320,000 realiza- 
tions. The error bars indicate one standard deviation, and 
the (red) horizontal line the asympotic value 3/2. 



system is half- filled (p = 1/2), and the hopping rates 
are equal in each direction transverse to the drive, but 
set totally asymmetric parallel to the drive. For a two- 
dimensional driven lattice gas (L x = L y = 100), we 
depict the effective exponent as defined in Eq. (fl~7| in 
Fig. II H and analogous results for the three-dimensional 
case are shown in Fig. 1121 Since we are at or above 
the critical dimension d c = 2, respectively, the asymp- 
totic value is £ = d/2. Indeed, the local exponent ( e ff(t) 
quickly relaxes to 1 for d = 2 and 3/2 for d = 3, and never 
surpasses these expected asymptotic values. The anoma- 
lously slow convergence and unusual finite-size behavior 
in one dimension are consequently caused by the strong 
out-of-cquilibrium correlations present only for d < d c . 



2.5 



1.5 



0.5 



d = 2 



~~m n 5T 



fj 



□ [] 



100 



200 300 
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FIG. 11. (Color online.) The effective exponent £ e ff(£), 
Eq. (|17p , for a two-dimensional driven lattice gas with L x — 
L y = 100 and p = 1/2. The data are averaged over 600,000 
realizations. The error bars indicate one standard deviation, 
and the (red) horizontal line the asympotic value 1. 



Up to this point, we have provided results only for 
driven lattice gases with zero backward hopping rate 
along the drive direction (TASEP in one dimension) . As- 
signing the forward hopping probability p and backward 
probability q, we can define the bias T = p — q = 2p — 1. 
As long as 1/2 < p < 1 and < T < 1 (ASEP), a macro- 
scopic particle current persists, the system remains out of 
equilibrium, and one expects qualitatively the same be- 
havior as for maximum bias or drive, p = T = 1. Thus, 
at long times, the temporal decay of the density auto- 
correlation function should be described by a power law 
with exponent £ = 2/3 in one dimension, for any non- 
zero value of the hopping bias T. However, note that this 
statement is valid at asympotically long times for suffi- 
ciently large systems, and a simple power-law decay will 
not necessarily be observable in finite systems with weak 
bias. Monte Carlo simulation data for one-dimensional 
driven lattice gases with various bias values are shown in 
Fig. Q21 As can be seen, the local exponent Ceff(i) relaxes 
increasingly slowly towards its asymptotic value 2/3 as 
the bias is reduced. The figure illustrates that for low 
bias and at short times, the density auto-correlation de- 
cay rather follow the symmetric exclusion process (SEP) 
relaxation, for which £ = 1/2. For example, for T = 0.2 
a distinct decrease in the local exponent is observed for 
t < 30, later followed by a very slow increase. 



IV. NON-EQUILIBRIUM RELAXATION 

A. Transitions between steady states 

Let us finally consider transient relaxation properties 
of driven lattice gases with exclusion. For example, start- 
ing from a steady state obtained with a specified set of 
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0.55 



25 50 75 100 125 150 175 200 
t (MCS) 

FIG. 13. (Color online.) Local auto-correlation exponent 
Ces{t) for one-dimensional driven lattice gases (ASEP) with 
different hopping biases T = 1.0 (TASEP), 0.6, 0.4, 0.2 (all 
ASEP), and 0.0 (vanishing drive, SEP); L = 4096 and p = 
1/2. The data are averaged over 5,000,000 realizations each. 
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FIG. 14. (Color online.) Finite-size scaling for the interface 
width of the one-dimensional growth model that corresponds 
to the TASEP, for system sizes L = 600, 800, and 1000, see 
Eq. (|27[) . The data are averaged over 50,000 realizations each. 



hopping probabilitcs p and q = 1 — p, we may wish to 
understand how the system responds to a sudden change 
of these probabilities to new values p' and q' = 1 — p' . 
However, recall that the non-equilibrium steady state 
probability distribution @ is independent of the hop- 
ping probabilities. This means that statistically the sta- 
tionary configurations before and after the sudden bias 
change are indistingishable. As a consequence, no slow 
relaxation processes follow the change of hopping proba- 
bilities. As we have confirmed with our Monte Carlo sim- 
ulations, any macroscopic observables such as the mean 
particle current assume their new stationary values es- 
sentially instantaneously after the bias reset. 



B. Strongly correlated initial conditions 

In order to construct initial conditions that generate 
non-stationary relaxation with broken time translation 
invariance, we consider the mapping of the TASEP to 
a surface growth problem in the KPZ universality class 
[4r| . In the growth model, one obtains slow relaxation to- 
wards the stationary state if the process is initiated with a 
flat surface h(x, 0) = everywhere, and roughening sub- 
sequently commences on small length scales ~ i 1 /*!! . In 
the TASEP picture, this initial state is accomplished at 
half- filling p = 1/2 by placing the particles alternatingly 
on every other site. The mean-square interface width of 
the corresponding finite one-dimensional growth model 
can be obtained from our TASEP Monte Carlo data on 
a ring of length L via Eq. ([1} , 



w(L,t) 2 = (h(x,t) 



!/=l 



x = l 



where w(L,0) 2 = 0. The standard finite-size scaling re- 
lation for the interface width reads [59 



w{L,t) = L a W{t/L z "), 



(27) 



(26) 



with a — 1/2 and zn — 3/2 in the KPZ universality class, 
which implies the initial growth law w(L,t) ~ t a l z \\ = 
t 1 / 3 , valid up to times t « L z w Our Monte Carlo 

simulation results nicely confirm the scaling relation (|27l) , 
as evidenced in Fig. [14] by the satisfactory data collapse 
for various system sizes with the expected scaling expo- 
nents (compare also the analogous mode-coupling results 
in Fig. 3 of Ref. 1 191 . and the simulation data shown 
in Fig. 1 of Ref. [57]). Our measurements yield an ini- 
tial growth exponent 0.30, slightly different from the ex- 
pected value 1/3, likely caused by sizeable corrections to 
scaling in our data. Our finite size scaling plot of the 
width matches the simulation data of Ref. [57| . 

Starting from the highly correlated initial configura- 
tion of alternatingly occupied sites in the TASEP that 
corresponds to a flat surface in the growth model rep- 
resentation, time translation invariance is broken during 
the slow approach to the stationary regime. The density- 
density or height-height correlation functions then de- 
pend explicitly on both time arguments t and t' . Hence- 
forth we shall refer to the second argument as waiting 
time s, and consider the short-time scaling limit s C t 
[32| . By means of a careful renormalization group analy- 
sis, Krech [3l[ established the following scaling form for 
the two-time height-height correlation function in this 
limit: 

Ci S (x,t,s«t) = ^ \x\ 2a F(t/x z «), (28) 

which generalizes Eq. (|22[) that is valid in the stationary 
regime. Here, the initial slip exponent 9 = (d + 4)/zii — 2 
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FIG. 15. (Color online.) Scaling plot for the density auto- 
correlation function for a one-dimensional driven lattice gas 
(TASEP) with length L = 1000 and density p = 1/2 in the 
non-equilibrium relaxation regime, see Eq. (|30|l . measured for 
waiting times s = 100, 160, and 220. The data are averaged 
over 60,000 realizations. 

can be expressed in terms of the standard dynamic scal- 
ing exponent z. This follows since the initial-time surface 
does not induce any novel singularities, as a consequence 
of the momentum dependence of the nonlinear vertex in 
the KPZ problem, and thus no additional rcnormaliza- 
tions are required. 

In one dimension, the renormalization group predicts 
6 = 4/3. For the TASEP density-density correlation 
function, we obtain by means of Eq. (|14[) 

S IS (x,t,s<^t) = (|)%| 2q - 2 /(*/^ii) 

= (f) 4/3 W- 1 /(*/« 3/3 )- (29) 

In order to arrive at the temporal behavior of the auto- 
correlation function, we must require that the scaling 
function /(£) ~ £~ 2 / 3 as its argument £ — ► oo, whence 

Sb(0, *,««*) = (f) 4/ V 2/3 - (30) 

Our Monte Carlo auto-correlation data displayed in 
Fig. 1151 obtained for various waiting times s, convinc- 
ingly confirm the scaling predicted by Eq. ([3"0|) for s<J. 

Finally, we explore the scaling properties of the two- 
time height and density auto-correlation functions in the 
physical aging regime. As shown for the KPZ growth 
model, no additional rcnormalizations, and hence no in- 
dependent new scaling exponents are required in the non- 
stationary relaxation regime in one dimension j3lj . Since 
this is essentially a consequence of the underlying con- 
served dynamics, we surmise that we can generalize the 
scaling laws (|15[1 and © to 

C(x,t,s) = b- 2a C(bx,b z H,b z "s), (31) 
S(x, t, s) = b d+A S(bx ± ,b 1+A x ll , b% b z s), (32) 



FIG. 16. (Color online.) Simple aging scaling plot, see 
Eq. (|19|) . for the density-density auto-correlation function 
for a one-dimensional driven lattice gas (TASEP) of length 
L = 1000, at density p = 1/2, obtained from Monte Carlo 
simulation data with different waiting times s = 100, 160, 
and 220. The data sets are averaged over 60,000 realizations 
each. 

with z = zn (1 + A) and 2z(l — a) = zu (d + A). Letting 
x — > and setting b = s -1 ^" , Eq. (pTTj) yields for the 
growth model height-height auto-correlation function 

C(0, t, s) = S 2a ' z * f c (t/s) = s 2/3 fc(t/s). (33) 

Recent Monte Carlo simulations have confirmed this sim- 
ple a ging scaling scenario [HH for the KPZ universality 
class [33l - l36| . In a similar vein, for the density-density 
auto-correlation function in the driven lattice gas, via 
matching b = s^ 1 ^ we obtain Eq. (fT9|) . with the scaling 
exponent £ given in Eq. (jlip . Figure [16] confirms the sim- 
ple aging scenario according to Eq. (|TT))) through the nice 
data collapse of our Monte Carlo simulation results for 
various waiting times in one dimension, where £ = 2/3. 
We have also performed driven lattice gas simulations in 
two and three dimensions with alternating particle ini- 
tial conditions on periodic lattices. As demonstrated in 
Figs. [T7] and [TH1 we obtain simple-scaling data collapse 
with the expected scaling exponent Q = 1 for d = 2 and 
C = 3/2for<i = 3, respectively. 

V. CONCLUSION 

Wc have performed extensive precise Monte Carlo sim- 
ulations of driven lattice gases with exclusion on periodic 
lattices, and studied in detail the behavior of the density 
correlations in one, two, and three dimensions. 

For the one-dimensional case, we have investigated in- 
teresting recurrence oscillations away from half-filling, 
and explained the minimum value in the density auto- 
correlation function observed for small systems. We also 
studied how varying the hopping bias in the ASEP af- 
fects the crossover to the long-time power law decay of 
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FIG. 17. (Color online.) Density auto-correlation function for 
a two-dimensional driven lattice gas of size L x — L y = 128 
for density p — 1/2 in the aging scaling regime for waiting 
times s — 50, 100, and 350. Each data set is averaged over 
160,000 realizations. 
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FIG. 18. (Color online.) Density auto-correlation function 
for a three-dimensional driven lattice gas of size L x — Ly — 
L z — 20 for density p = 1/2 in the aging scaling regime for 
waiting times s — 10, 20, and 30. Each data set is averaged 
over 4,200,000 realizations. 



the density auto-correlation function. While our large- 
scale simulations confirm standard finite-size scaling be- 
havior, the effective exponent that describes the tempo- 
ral decay of the auto-correlations approaches the asymp- 
totic value exceedingly slowly, even for the TASEP, and 
displays very unusual finite-size corrections to scaling. 
Our Monte Carlo simulation results confirm all expected 
non-trivial scaling exponents for the driven lattice gas, in 
accord with the scaling properties of equivalent growth 
models in the KPZ universality class. This includes pre- 
dictions for the simple aging scaling scenario during the 
non-stationary relaxation from strongly correlated initial 
conditions. Our data thus support the assertion that 
in driven lattice gas with exclusion, the entire universal 
scaling regime is governed by a single non-trivial expo- 
nent A (or equivalcntly, z» or £). Specifically, no novel 
scaling exponent is required to capture non-equilibrium 
aging properties in this system. We have moreover nu- 
merically determined several different universal scaling 
functions. 

We extended our Monte Carlo simulations to two and 
three dimensions. Aside from logarithmic corrections at 
the critical dimension d c = 2, we have confirmed mean- 
field scaling behavior for the density auto-correlation 
function in both the stationary long-time and the non- 
stationary aging regimes. 
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